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Abstract 

The contact dynamics method (CD) is an efficient simulation technique of dense 
granular media where unilateral and frictional contact problems for a large number 
of rigid bodies have to be solved. In this paper we present a modified version of the 
contact dynamics to generate homogeneous random packings of rigid grains. CD 
is coupled to an external pressure bath, which allows the variation of the size of a 
periodically repeated cell. We follow the concept of the Andersen dynamics and show 
how it can be applied within the framework of the contact dynamics method. The 
main challenge here is to handle the interparticle interactions properly, which are 
based on constraint forces in CD. We implement the proposed algorithm, perform 
test simulations and investigate the properties of the final packings. 
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1 Introduction 



Computer simulation methods have been widely employed in recent years to 
study the behavior of granular materials. Among the numerical techniques, 
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discrete element methods, including soft particle molecular dynamics (MD) 
[1,2], event-driven (ED) [3,4] and contact dynamics (CD) [5,6,7,8], constitute 
an important class where the material is simulated on the level of particles. 
In such algorithms the trajectory of each particle is calculated as a result 
of interaction with other particles, confining boundaries and external fields. 
The differences between the discrete element methods stem from the way how 
interactions between the particles are treated, which leads also to different 
ranges of applicability. 

In low density granular systems, where interactions are mainly binary col- 
lisions, the event-driven method is an efficient technique. The particles are 
modeled as perfectly rigid and the contact duration is supposed to be zero. 
The handling of dense granular systems, where the frequency of collisions is 
large or long-lasting contacts appear, becomes problematic in ED simulations 
[9,10]. 

In case of dense granular media the approach of soft particle molecular dy- 
namics is more favorable and widely used. In MD, the time step is usually 
fixed and the original undeformed shapes of the particles may overlap during 
the dynamics. These overlaps are interpreted as elastic deformations which 
generate repulsive restoring forces between the particles. Based on this inter- 
action, which is defined in the form of a visco-elastic force law, the stiffness of 
the particles can be controlled. When the stiffness is increased MD simulations 
become slower since the time step has to be chosen small enough so that the 
velocities and positions vary as smooth functions of time. 

The contact dynamics method considers the grains as perfectly rigid. Therefore 
no overlaps between the particles are expected and they interact with each 
other only at contact points. The contact forces in CD do not stem from 
visco-elastic force laws but are calculated in terms of constraint conditions 
(for more details see Sec. 2). This method has shown its efficiency in the 
simulation of dense frictional systems of hard particles. 

Packings of hard particles interacting with repulsive contact forces are ex- 
tensively used as models of various complex many-body systems, e.g. dense 
granular materials [11], glasses [12], liquids [13] and other random media [14]. 
Jamming in hard-particle packings of granular materials has been the subject 
of considerable interest recently [15,16]. Furthermore hard-particle packings, 
and especially hard-sphere packings, have inspired mathematicians and been 
the source of numerous challenging theoretical problems [17], from which many 
are still open. 

Real systems in the laboratory and in nature contain far too large number 
of particles to model the whole system in computer simulations. Due to the 
limited computer capacity the simulations are often restricted to test a small 
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mesoscopic part of a large system. Typically, the studies are focused to a 
local homogeneous small piece of the material inside the bulk far from the 
border of the system. Therefore simulation methods are required that are 
able to generate and handle packings of hard particles without side effects of 
confining walls. 

The usual simulation methods of dense systems involve confining boxes where 
the material is compactified by moving pistons or gravity. However, the prop- 
erties of the material differ in the vicinity of walls and corners of the confining 
cell from those in the bulk far from the walls. The application of walls in 
computer simulations leads to inhomogeneous systems due to undesired side 
effects (e.g. layering effect). Moreover, the structure of the packings becomes 
strongly anisotropic in these cases due to the orientation of walls and special 
direction of the compaction. For studies where such anisotropy is unwanted 
other type of compaction methods are needed. 

In this paper, we present a compaction method where boundary effects are 
avoided due to exclusion of side walls. This simulation method is based on the 
contact dynamics algorithm where we applied the concept of the Andersen 
dynamics [18], which enables us to produce homogeneous granular packings of 
hard particles with desired internal pressure. The compaction method involves 
variable volume of the simulation cell with periodic boundary conditions in 
all directions. 

This paper is organized as follows. First we present some basic features of CD 
method in Section 2. Then, Section 3 describes the equations of motion for a 
system of particles with variable volume. In Section 4 we present a modified 
version of CD with coupling to an external pressure bath. In Section 5 we 
report the results of some test simulations. Section 6 concludes the paper. 



2 Non smooth contact dynamics 

Contact dynamics (CD), developed by M. Jean and J. J. Moreau [6,7,8], is a 
discrete element method in the sense that the time evolution of the system is 

— * 

treated on the level of individual particles. Once the total force Fj and torque 

— * 

Ti acting on the particle i is known, the problem is reduced to the integration 
of Newton's equations of motion which can be solved by numerical methods. 
Here we use the implicit first order Euler scheme: 

1 

v i (t + At)=v i (t) + —F i (t + At)At (1) 
mi 

fi(t + At) = n(t) + Vi(t + At) At , (2) 
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which gives the change in the position r*j and velocity Vi of the center of mass 
of the particle with mass mj after the time step At. At is chosen so that 
the relative displacement of adjacent particles during one time step is small 
compared to the particle size and to the radius of curvature of the contacting 
surfaces. Corresponding equations are used also for the rotational degrees 
of freedom, describing the time evolution of the orientation and the angular 
velocity uji caused by the new total torque Tj(t + At) acting on the particle i. 

The interesting part of the CD method is how the interaction between the 
particles are handled. For simplicity we assume that the particles are nonco- 
hesive and dry, we exclude electrostatic and magnetic forces between them 
and consider only interactions via contact forces. The particles are regarded 
as perfectly rigid and the contact forces are calculated in terms of constraint 
conditions. Such constaints are the impenetrability and the no-slip condition, 
i.e. the contact force has to prevent the overlapping of the particles and the 
sliding of the contact surfaces. This latter condition is valid only below the 
Coulomb limit of static friction, which states that the tangential component 
of a contact force R can not exceed the normal component times the friction 
coefficient /i: 



If the friction is not strong enough to ensure the no-slip condition the contact 
will be sliding and the tangential component of the contact force is given by 
the expression 



where v$~ stands for the tangential component of the relative velocity between 
the contacting surfaces. In the CD method the constraint conditions are im- 
posed on the new configuration at time t + At, i.e., the unknown contact forces 
are calculated in a way that the constraints conditions are fulfilled in the new 
configuration [8]. This is the reason why an implicit time stepping is used. 

In order to let the system evolve one step from time t to t + At one has 
to determine the total force and torque acting on each particle which may 
consist of external forces (like gravity) and contact forces from neighboring 
particles. Let us suppose that all the unknown contact forces are already 
determined except for one force between a pair of particles already in contact 
or with a small gap between them. Here we explain briefly how the constraint 
conditions help to determine the interection between these two particles. A 
detailed description of the method can be found in [8] . 

The algorithm starts with the assumption that the contact force we are search- 
ing for is zero and checks whether this leads to an overlap of the undeformed 
shapes of the two particles after one time step At. This is done based on the 
time stepping [Eq. (1)]: The external forces and other contact forces provide 



(3) 



Rt — —[i'R, 
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(4) 
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Fi(t + At) and Tj(t + At) for both particles thus the new relative velocity of 
the contacting surfaces rel ' free can be calculated. Here we use the term con- 
tacting surfaces for simplicity thought the two particles are not necessarily in 
contact. There might be a positive gap g between them, which is the length 
of the shortest line connecting the surfaces of the two particles (Fig. 1). We 
will refer to the relative velocity of the endpoints of the line as the relative 
velocity of the contact and denote the direction of the line by the unit vector 
n. In the limit of a real contact g is zero and n becomes the contact normal. 
Negative gap has the meaning of an overlap. The superscript free in t? rel > free 
denotes that the relative velocity has been calculated assuming no interaction 
between the two particles. We use the sign convention that negative normal 
velocity (n ■ ■i? rel < 0) means approaching particles. 

The new value of the gap (after one time step) is estimated by the algorithm 
based on the current gap g and the new relative velocity ■y ¥cl,frce according to 
the implicit time stepping. If the new gap is positive: 

g + {f el ' free • nAt > (5) 

then the zero contact force (no interaction) is accepted because no contact 
is formed between the two particles. However, if the estimated new gap is 
negative then a contact force has to be applied in order to avoid the viola- 
tion of the constraint conditions. Generally, one expects the following relation 
between the unknown new contact force R and the unknown new relative 
velocity {F cl : 

^ = ^M(if el ' free -{f el ), (6) 
where M is the mass matrix that describes the inertia of the contact, i.e. M i? 

— # 

is the relative acceleration of the contacting surfaces due to the contact force R. 
The mass matrix M depends on the shape, mass and moment of inertia of the 
two particles. On one hand, the interpenetration of the two rigid particles has 
to be avoided, which gives the following constraint for the normal component 
of {F el : 

g + v" cl ■ nAt = 0. (7) 

On the other hand, the tangential component of -iF cl has to be zero in order 
to ensure the no-slip condition 

^ = 0. (8) 
The required contact force that fulfills Eqs. (7) and (8) then reads 

R = ^M(^-n + {f el ' free ). (9) 
At y At ' w 

This contact force is acceptable only if it fulfills the Coulomb condition [Eq. (3)]. 
Otherwise we can not exploit the non-slip contact assumption. In this case, 
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vl el is not zero, the contact slides and the contact force has to be recalculated. 
Eqs. (6) and (7) then provide 

R = ^-M(-?-n - + {f el > free ), (10) 
At v At * 1 v 1 

where the number of unknowns (components of R and vl° l ) exceeds the number 

— # 

of equations. In order to determine the contact force R one has to solve Eq. (10) 
together with Eq. (4). 

It is recommended to use g pos = max(g,0) instead of g in Eqs. (9) and (10). 
The gap size g should always be non-negative and using g pos apparently makes 
no difference [7]. However, due to the inaccuracy of the calculations small 
overlaps can be created between neighboring particles. If g instead of g pos is 
used then these overlaps are eliminated in the next time step by imposing 
larger repulsive contact forces to satisfy Eq. (7), which pumps kinetic energy 
into the system. Using g pos instead of g eliminates this artifact on the cost 
that an already existing overlap is not removed (which then serves to check 
the inaccuracies of the simulation [19]), only its further growth is prevented. 
Regarding the above mentioned points, we rewrite the equations (9) and (10) 
as 

_ 1 p° s 

R = — M(V- n + ^ el,free ) and (11) 
At v At 1 v ' 



R = M(^—n - vf + £* el > free ). (12) 
At v At * 1 v ' 

A flowchart of the single contact force calculation is given in Fig. 2. So far 
we have explained only how the CD algorithm determines a single existing or 
incipient contact, based on the assumption that all the surrounding contact 
forces are known. However, in a dense granular media, many particles contact 
simultaneously and form a contact network. In this case, a contact force cannot 
be evaluated locally, since it depends on the adjacent contact forces which are 
also unknown. To find a globally consistent force network at each time step, 
an iterative scheme is applied in CD. 

At each iteration step, all contacts are chosen one by one and the force at the 
contact is updated according to the scheme shown in Fig. 2. The update is 
sequential, i.e., the freshly updated contact force is stored immediately as the 
current force and then a new contact is chosen for the next update. 

After one iteration step, constraint conditions are not necessarily fulfilled for 
each contact. In order to find a global solution the iteration process has to 
be repeated several times until the resulting force network converges. The 
convergence of the iteration process is smooth, i.e., the precision of the solution 
increases with the number of iterations Nj. Higher Nj provides more precise 
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solution but also requires more computational effort. The CD method can 
be used with a constant number of iterations in subsequent time steps [19,8] 
or with a convergence criterion that prescribes a given precision to the force 
calculation [6,7,8]. In this latter case, the number of iterations AT/ varies from 
time step to time step. 

After the new force network is determined with a prescribed precision, the 
system evolves at the end of the time step according to the time-stepping 
scheme described at the beginning of this section. It is important to note that 
choosing small Ni and/or large time step causes systematic errors of the force 
calculation which lead to a spurious soft particle behavior [19] in spite of the 
original assumption of perfect rigidity. 

To conclude this section, we present briefly the scheme of the solver. 

t := t + At (time step) 
Evaluating the gap g for all contacts 
Nj := Nj + 1 (iteration) 
k := k + 1 (contact index) 

Evaluating t^ el ' free then (according to the flowchart in Fig. 2) 
Convergence test for contact forces 
Time-stepping for velocities and positions of all particles (using Eqs. (1) and (2)) 

3 The equations of motion at constant external pressure 

In the simulation of granular materials, it is often desirable to investigate 
systems which are not surrounded by walls and to apply periodic boundary 
conditions in all directions. It is a nice feature of periodic boundary condi- 
tions that they make points of the space equivalent, the boundary effects are 
eliminated. That way the bulk properties of the material can be studied more 
easily. However, the application of an external pressure becomes problematic 
since the total volume is fixed and the system cannot be compressed by pistons 
or moving walls. 

In order to overcome this problem, but at the same time keep the advantageous 
periodic boundaries, Andersen [18] proposed a method for molecular dynamics 
simulations. Here, we recall his method briefly as its basic ideas will be used 
later on in this paper. 

According to the Andersen method the boundaries are still periodically con- 
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nected in all directions and no walls are present, but the volume of the system 
is a dynamical variable which evolves in time due to interaction with an ex- 
ternal pressure bath. When a system of iV atoms is compressed or expanded 
it is done in an isotropic and homogeneous way: The distances between the 
atoms are rescaled by the same factor regardless of the relative or absolute 
positions. 

Let us give the equations of motion of a system with particle positions fi, f 2 , fjv 
in a cubic volume V (each component of fj is between and V 1 ^): 

dm _ p^) i dhiV(t) 

m =m .^jEm, (14) 
M^m =Mt) . P „ t=AP(t) . (15) 

Eq. (13) describes the change in the position r*;. The first term on the right is 
the usual one, the momentum pi divided by the mass of the ith particle. The 
last term is the extension by the Andersen method that rescales the position 
according to the relative volume change. 

Eq. (14) provides the time evolution of the momentum due to two terms. The 

— * 

first one is the usual total force Fi acting on the ith particle which originates 
from the interaction with other particles and/or from external fields. The 
additional term leads to further acceleration of the particle if the volume is 
changing. E.g., if the system is compressed the kinetic energy of the particles 
is increased due to the work done by the compression. The energy input is 
achieved by rescaling all particle momenta regardless of their positions. This 
is in contrast to usual pistons where the energy enters at the boundary. 

Eq. (15) can be interpreted as Newton's second law that governs the change 
of the volume. It describes the time evolution of an imaginary piston which 
has the inertia parameter M v and is driven by the generalized force AP(£) = 
Pm{t) — -fext- This latter is the pressure difference between the external pressure 
P cxt imposed by the pressure bath and the internal pressure of the system 
P m (t). The pressure difference AP(t) drives the system towards the external 
pressure, the sensitivity of the system to this driving force is controlled by the 
inertia parameter M v . 

In the limit of infinite inertia M v — > oo together with the initial condition 
dV(t )/dt = the volume of the system remains constant and Eqs. (13) 
and (14) correspond to the usual Newtonian dynamics of the particles. 

In order to get more insight into the Andersen dynamics let us consider a 
simple example of a system of non-interacting particles with all Fi(t) = 0. Ini- 
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tially, the velocities and the volume velocity dV(to)/dt are set to zero. Because 
the internal pressure P ni is zero the system with finite inertia M v and under 
external pressure P ext > will start contracting according to Eq. (15). The 
acceleration of the particles [Eq. (14)] remains zero during the time evolution; 
One might say that the particles are standing there all the time. However, 
the distances between them are decreasing because of the contraction of the 
"world" around them. This is caused by the second term on the right hand 
side of Eq. (13) while the first term remains zero. 

This suggests the picture of an imaginary background membrane that con- 
tracts or dilates homogeneously together with the volume and carries the par- 
ticles along. The velocity of this background at position r*; is given by A(i)fj 
where A is the dilation rate defined by the rate of the relative change in the 
system size L: 

L(t) 1 dlnVjt) 
A(t) = W) = D^r~ (16) 
and D is the dimension of the system. Then the right hand side of Eq. (13) can 
be interpreted as the sum of two velocities: the second one is the velocity of 
the background at the position of the particle and the first one is the intrinsic 
velocity of the particle measured compared to the background. The sum of 
these two forces gives the changing rate of the absolute position r*j. In the 
rest of the paper the velocity Vj will refer always to the intrinsic velocity. We 
rewrite Eq. (13) in the following form 

d ^ = v t (t) + mm. (17) 



Next we turn to the modelling of granular systems. Our goal is to achieve static 
granular packings that are compressed from a loose gas-like state. Here again 
it is advantageous to exclude confining walls and in order to apply pressure 
and achieve contraction of the volume we will use the concept of the Andersen 
method. However, the equations of motion will be slightly changed in order 
to make them suit better to our goals. 

In granular materials the interaction between the particles are dissipative. 
When the material is poured into a container or is compressed by a piston the 
particles gain kinetic energy due to the work done by gravity or the piston. All 
this energy has to be dissipated (turned into heat) by the interactions between 
the particles before the material can settle into a static dense packing of the 
particles. This relaxation process requires a massive computational effort when 
large packings are modeled in computer simulations. One encounters the same 
problem if the Andersen dynamics is applied straight to granular systems. The 
relaxation time can be reduced if the second term on the right hand side of 
Eq. (14) is omitted, because then the total amount of energy pumped into the 
system is reduced. In this case the particles are accelerated only by the forces 
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Fi but they receive no additional energy due to the decreasing volume. Thus 
the following equation will be applied for granular compaction here 



which results in a more effective relaxation rather than Eq. (14). This change 
is advantageous also from the point of view of momentum conservation. If 
the system is compactified by using Eq. (14) from an initial condition where 
the total momentum of the particles is non-zero, then this momentum will be 
increased inverse proportionally to the size of the system L. Thus the total 
momentum e.g. due to initial random fluctuations is magnified which can lead 
to non-negligible overall rigid body motion of the final static packing. This is 
in contrast to Eq. (18) which provides momentum conservation in the absence 
of external fields. 



Concerning the equation that describes the time evolution of the system size 
we find it more convenient to control A instead of dV/dt. This is actually not 
an important change and leads to very similar dynamics. Our third equation 
reads 

M x ^± = AP(t). (19) 



The equations of motion (17), (18) and (19) describe an effective compaction 
dynamics for granular systems, they are able to provide static packings under 
the desired pressure P cxt and if they are restricted to the limit of M\ — > oo 
we receive back the classical Newtonian dynamics. 

Of course, in order to close the equations we need to define interactions be- 
tween the particles. The interparticle forces provide in Eq. (18) and they 
are also needed to evaluate the inner pressure P m . 

The stress tensor a a p is not apriori spherical in granular materials. The average 
a a p of the system is determined by the interparticle forces [20] and the particle 
velocities as 

I N c N 

Va(3 = 77 F k,Jk,(3 + m i v i,a v i,P&ap) , (20) 
V k=l i=l 

where N and N c denote the number of particles and the number of contacts, 
respectively. If two contacting particles at contact k are labelled by 1 and 2, 

— # — * 

then Fk is the force exerted on particle 2 by particle 1, and the vector 4 is 
pointing from the center of mass of particle 1 to that of particle 2 where peri- 
odic boundary conditions and nearest image neighbors are taken into account. 
Thus Ik is the minimum distance between particles 1 and 2: 

l k = \f k \ =m i n \f 2 -f 1 + V 1/3 a\, (21) 
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where a is an integer-component translation vector. The inner pressure is then 
given by the trace of the stress tensor divided by the dimension of the system 

1 n c N 

P in = 7777 E ' h + E * 3] > ( 22 ) 

uv fc=l 1=1 
which has the meaning of an average normal stress. 

The implementation of the above method in computer simulations is straight- 
forward if the interparticle forces are functions of the positions and velocities 
of the particles, e.g., in soft particle MD simulations. The implementation is 
less trivial for the case of the contact dynamics method where interparticle 
forces are constraint forces. We devote the next section to this problem. 



4 Contact dynamics with coupling to an external pressure bath 



In this section we present a modified version of the contact dynamics algorithm 
which enables us to perform CD simulations when the system is in contact 
with an external pressure bath. According to Sec. 3 let us suppose that the 
system is subjected to a constant external pressure P ext and its time evolution 
is given by Eqs. (17), (18) and (19). 

Here we will follow the description of the CD method given in Sec. 2 and 
discuss the required modifications. Once the force calculation process is com- 
pleted, the implicit Euler integration can proceed one time step further. Now 
the time stepping has to involve also the equations of motion of the system 
size. By discretizing the Eqs. (16) and (19) in the same implicit manner as for 
the particles [Eqs. (1) and (2)] we obtain the new values of the system size L 
and the dilation rate A: 

X(t + At) = X(t) + AP( ^+ At) At, (23) 

L{t + At) = L{t) [1 + X(t + At) At] (24) 

where the "velocity" X(t) and the "position" L(t) are updated by the new 
"force" AP(t + At) and by the new "velocity" X(t + At), respectively. 

The discretized equations governing the translational degrees of freedom of 
the particles [Eqs. (1) and (2)] are rewritten according to Eqs. (17) and (18) 
in the following form: 

Vi(t + At) = Vi(t) + — Fi(t + At) At and (25) 

fi(t + At) = fi(t)[l + X(t + At) At] + vft + At) At. (26) 
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The time stepping for the rotational degrees of freedom remains unchanged 
because the dilation (contraction) of the system has no direct effect on the 
rotation of the particles. 

In the CD method, as we explained in Sec. 2 the particles are perfectly rigid 
and are interacting with constraint forces, i.e. those forces are chosen between 
contacting particles that are needed to fulfill the constraint conditions. E.g. the 
contact force has to prevent the interpenetration of the contacting surfaces. 

If an external pressure bath is used then the calculation of the constraint 
forces has to be reconsidered because the relative velocity of the contacting 
surfaces is influenced by the variable volume. When the system is dilating 
or contracting, particles gain additional relative velocities compared to each 

— * — * 

other. For a pair of particles, this velocity is XI where I is the vector connecting 
the two centers of mass. The same change appears in the relative velocity of 
the contacting surfaces as the size of the particles is kept fixed. If this change 
led to interpenetration then it has to be compensated by a larger contact force. 
It may also happen that existing contacts open up due to expansion of the 
system resulting in zero interaction force for those pair of particles. 

— * 

In the calculation of a single contact force, the relative velocity XI (i.e. the 
contribution of the changing system size) has to be added to ■y' rel ' free . The 
new relative velocity of the contact assuming no interaction between the two 
particles ■y' rel ' free is calculated here in the same way as in Sec. 2, i.e., based on the 
intrinsic velocities of the particles. Thus the effect of the dilation/contraction 
of the system is not taken into account in ^ el > free . Therefore one has to replace 
{jrei,free with (■y Tel ' free + XI) in all equations of Sec. 2 in order to impose the 
constraint conditions properly. Let us first suppose that the system has infinite 
inertia (Ma = oo) thus the dilation rate A is constant. In this case the modified 
equations of the force update (containing already the term •u rel ' free +AZ) provide 
the right constraint forces at the end of the iteration process. These forces will 
alter the relative velocity (^ el ' free + XI) in such a way that the prescribed 
constrain conditions will be fulfilled in the new configuration at t + At. 

More consideration is needed if finite inertia M\ is used and the dilation rate 
A is time-dependent. The problem is that in order to calculate the proper 
contact force one has to know the new dilation rate. The new dilation rate, 
however, depends on the new value of the inner pressure [Eq. (23)] which, in 
turn, depends on the new value of the contact forces. This problem can be 
solved by incorporating A and P m into the iteration process. Instead of using 
the old values X(t) and P m (t) during the iteration, we always use the expected 
values A* and P* n . These represent our best guess for the new dilation rate 
X(t + At) and for the new inner pressure P m {t + At). P* n is defined based on 
the current values of the contact forces F£ during the force iteration. Whenever 
a contact force is updated we recalculate the expected inner pressure P* n . With 
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the help of the current contact forces F£ we can determine the total forces 
acting on the particles and then, following Eq. (25) we obtain also the expected 
new velocities of the particles v*. The expected inner pressure, according to 
Eq. (22), then reads: 

i n c N 

K = ™7E ■ h + • (27) 

1JV fc=l 1=1 

Of course, there is no need to recalculate all the terms in Eq. (27) in order to 
update P* n . When the force at a single contact is changed it affects only three 
terms: one due to the force itself and two due to the velocities of the contacting 
particles. In order to save computational time, only the differences in these 
three terms have to be taken into account when P* n is updated. Following 
Eq. (23), we obtain also the corresponding value of the expected dilation rate: 

A* = A(t) + P ' n ~ Pcxt At (28) 
M x 

This way, A* and P* n are updated many times between two consecutive time 
steps (in fact they are updated NjN c times) but in turn A* and P* n are always 
consistent with the current system of the contact forces. At the end of the 
iteration process P* n and A* provide not just an approximation of the new 
inner pressure and new dilation rate but they are equal to P m {t + At) and 
A(t + At), respectively. 

To complete the algorithm, we list here also the equations that are used for 
the force calculation of a single contact. The inequality (5) is replaced by 

g + (^ cl > frcc + A* T) • nAt > 0, (29) 

i.e. there is no interaction between the two particles if the inequality is satisfied. 
Otherwise we need a contact force. The force, previously given by Eq. (11), 
that is required by a sticking contact is 

R = M(^—n + A* / + {F el < free ). (30) 
At v At ' v ' 

This force again has to be recalculated according to a sliding contact if R in 
Eq. (30) violates the Coulomb condition: 

R = M(^—n - v[ cX + A* 1 + {f el ' free ), (31) 
At v At * 1 v 1 

which replaces the original equation (12). 

Except the above changes, the CD algorithm remains the same. In each time 
step the same iteration process is applied in order to reach convergence of the 
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contact forces. After the iteration process we apply Eqs. (23)- (26) to complete 
the time step. 

The scheme of the solver for the modified version of CD can be presented as 

t = t + At (time step) 
Evaluating the gap g for all contacts 
Nj = Nj + 1 (iteration) 
k = k + 1 (contact index) 
Evaluating v k e1, ree 

Evaluating R k (using Eqs. (30) and (31)) 

Evaluating P* n (Eq. (27)) then A* (Eq.( 28)) 
Convergence test for contact forces 
Time-stepping for the dilation rate and the system size (using Eqs. (23) and (24)) 
Time-stepping for velocities and positions of all particles (using Eqs. (25) and (26)) 

In the next section we will present some simulations with the above method. 
We will test the algorithm and analyse the properties of the resulting packings. 

As an alternative to this fully implicit method we considered another possibil- 
ity to discretize the Eqs. (16)-(19) in the spirit of the contact dynamics and, 
at the same time, impose the constraint conditions on the new configuration. 
The main difference is that the new value of the inner pressure P«„(t + At) is 
determined based on the old velocities Vi(t) and not on the new ones Vi(t+At), 
while the contribution of the forces are taken into account in the same way, i.e. 
the new contact forces F k (t + At) are used in Eq. (22). Therefore this version 
of the method is only partially implicit, however, the constraint conditions and 
the force calculation [Eqs. (29)-(31)] can be applied in the same way. Only, 
the expected values A* and P* n has to be changed consistently with the new 
pressure Pi n {t + At): 

i N c N 

K = 7ST7 E P k ■ ^ + E rrum ' (32) 
uv k=i i=i 

and then this modified P* n is used to determine the expected dilation rate A* 
with the help of Eq. (28). Again here, P* n and A* are calculated anew after 
each force update during the iteration process and their last values equal the 
new pressure Pi n (t + At) and the new dilation rate X(t + At). 

We implemented and tested this second version of the method and found 
that the constraint conditions are handled here also with the same level of 
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accuracy. Although the second method is perhaps less transparent than the 
fully implicit version, for practical applications it seems to be more useful. 
First, the second version of the method is easier to implement into a program 
code, second, it turned out to be faster by 25% in our test simulations. The 
improvement of the computational speed originates from the smaller number 
of the operations. One does not have to handle the expected particle velocities 
v* and the recalculation of P* n is more simple as the change of a contact force 
F£ affects only one term in the Eq. (32). 



5 Numerical results 

We perform numerical simulations using the CD algorithm with the fully im- 
plicit pressure bath scheme of section 4. This algorithm has been used to study 
mechanical properties of granular packings in response to local perturbations 
[21,22]. Here, the main goal is to show that the algorithm works indeed in 
practical applications and to test the method from several aspects; We in- 
vestigate how the simulation parameters influence the required CPU time and 
the accuracy of the simulation. Such parameters are the external pressure P ext , 
the inertia parameter M\ and the computational parameters, like the number 
of iterations per time step Nj and the length of the time step At. We also 
analyse the properties of the resulting packings. 

Here, we report only simulations of two-dimensional systems of disks, where 
the behavior is very similar to that we found for spherical particles in three- 
dimensional systems. Length parameters, the time and the two-dimensional 
mass density of the particles are measured in arbitrary units of Iq, to and po, 
respectively. The samples are polydisperse and the disk radii are distributed 
uniformly between 0.8 and 1.2, thus the average grain radius is 1. The material 
of the grains has unit density and the masses of the disks are proportional to 
their areas. In this section we have one reference system that contains 100 
disks. The interparticle friction coefficient is set to 0.5. The value of other 
parameters are: iVj = 100, At = 0.01, P cxt = 1 (this latter is expressed in 
units of p lo 2 /t 2 ) and the inertia M\ = 100 (in units of po^o 2 )- Throughout 
this section, we either use these reference parameters or the modified values 
will be given explicitly. Usually, we will vary only one parameter to check its 
effect while other parameters are kept fixed at their reference values. 

In each simulation, we start with a dilute sample of nonoverlapping rigid 
disks randomly distributed in a two dimensional square-shaped cell. No con- 
fining walls are used according to the boundary conditions specified in Sec. 4. 
Gravity and the initial dilation rate are set to zero. Due to coupling to the 
external pressure bath the dilute system starts shrinking. As the size of the 
cell decreases, particles collide, dissipate energy and after a while a contact 
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force network is formed between touching particles in order to avoid interpen- 
etrations. The contact forces build up the inner pressure P- m which inhibits 
further contraction of the system. Finally, a static configuration is reached in 
which P in equals P cxt and mechanical equilibrium is provided for each parti- 
cle. Technically, we finish the simulation when the system is close enough to 
the equilibrium state: We apply a convergence threshold for the mean veloc- 
ity fmcan and mean acceleration a mcan of the particles (which are measured 
in units lo/t and /o/to 2 , respectively). Only if both f mea n and a mean become 
smaller than the threshold 10~ 10 we regard the system as relaxed and stop the 
simulation. 

The typical time evolution can be seen in Fig. 3 where we show the compaction 
process in the case of the reference system. Fig. 3 (left) implies that the mag- 
nitude of A grows linearly in the beginning when the inner pressure is close to 
zero. The negative value of the dilation rate indicates contraction which be- 
comes slower after the particles build up the inner pressure [Fig. 3(middle)]. 
The fluctuations in P m are due to collisions of the particles. In the final stage 
of the compression A goes to zero, P in converges to the external pressure and 
the size of the system reaches its final value [Fig. 3(right)]. 

Next we investigate how the required CPU time of the simulation is affected 
by the various parameters. All simulations are performed with a processor 
Intel(R) Core(TM)2 CPU T7200 @ 2.00GHz and the CPU time is measured 
in seconds. Figure 4 reveals that the variation of P cx t, At, Nj and M\ have 
direct influence on the required CPU time. The final packing is achieved with 
less computational expenses if larger P ext , larger At or smaller Nj is used. The 
role of system inertia M\ is more complicated. M\ reflects the sensitivity of 
the system to the pressure difference p n — P ext . If the level of the sensitivity is 
too small or too large, the simulation becomes inefficient. It is advantageous to 
choose the inertia M\ near to its optimal value which depends on the specific 
system (e.g. on the number and mass of the particles) [23]. 

Regarding the efficiency of the computer simulation, not only the computa- 
tional expenses play an important role but the accuracy of the simulation is 
also essential. Here we use the overlaps of the particles as a measure of the 
inaccuracy of the simulation (see Sec. 2). In an ideal case there would be 
no overlaps between perfectly rigid particles. Fig. 5 shows the mean overlaps 
measured in the final packings. It can be seen for the parameters P ex t, At and 
Nj that the reduction of the computational expenses at the same time leads 
also to the reduction of the accuracy of the simulation. In Fig. 6 (left) we plot 
the CPU time versus the mean overlap for these three parameters. The data 
points collapse approximately on the same master curve which tells us that 
the computational expenses are determined basically by the desired accuracy; 
smaller errors require more computations. The efficiency of the simulation is 
approximately independent of the parameters P ext , At and Nj. The situation 
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is, however, different for the inertia of the system M\ . First, it has relatively 
small effect on the accuracy of the simulation (Fig. 5). Variation of M\ by 7 
orders of magnitude could hardly change the mean overlap of the particles. 
Second, M A affects strongly the efficiency of the simulation which is shown 
clearly by Fig. 6(right). If M\ is varied then larger computational expense is 
not necessarily accompanied with smaller errors. In fact, in the whole range 
of M\ studied here, the fastest simulation turned out to be the most accurate 
one [open circle in Fig. 6(right)]. 

Next we turn to the question whether the parameters of the simulation used in 
the compaction process influence the physical properties of the final packing. 
There are many ways to characterize static packings of disks. Here we test only 
one quantity, the frequently used volume fraction. The volume fraction gives 
the ratio between the total volume (total area in 2D) of grains and the volume 
(area) of the system. Fig. 7 shows the volume fraction of the same packings 
that were studied already in Fig. 5. It can be seen that the volume fraction 
remains approximately unchanged under the variation of the four parameters 
-Pext, Nj, At and M\ . This is except for one data point for large time step 
At, where the simulation is very inaccurate. The corresponding mean overlap 
(Fig. 5) is comparable to the typical size of the particles which is the reason 
why the volume fraction appears to be much larger. 

Finally, we investigate the inner structure of the resulting random packings. 
For that, we study larger samples with N = 1000 particles, otherwise, the 
default parameters are used during the compaction. To suppress random fluc- 
tuations, we produce 5 different systems and all quantities reported hereafter 
represent average values over these systems. In Fig. 8 (left) we study the angu- 
lar distribution of the contact normals and find that it is very close to uniform. 
However, there is a small but definite deviation (around 3%): the density of 
the contact normals are slightly larger parallel to the periodic boundaries. In 
this sense the packing is not completely isotropic. Although the effect is very 
small, the orientation of the boundaries can be observed also in such local 
quantities like the direction of the contacts. 

In connection to the question of the isotropy we checked also the global stress 
tensor a. In the original frame a reads 

I 1.00909 -0.01334 \ 

(33) 

^-0.01332 0.99091 j 

This stress is isotropic with good approximation. The diagonal entries are 
close to 1 which equals P ex t while the off-diagonal elements are approximately 
zero. Compared to the unit matrix, the elements deviate around 1% of P e xt- 

The final packings are expected to be homogeneous as all points of the space 
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are handled equivalently by the compaction method. Apart from random fluc- 
tuations we do not observe any inhomogeneity in our test systems. As an 
example, we show the spatial distribution of the contacts in Fig. 8 (right), 
where the density is approximately constant. 



6 Concluding remarks 

In this work we have proposed and tested a simulation method to produce 
homogeneous random packings in the absence of confining walls. We com- 
bined the contact dynamics algorithm with a modified version of the Ander- 
sen method to handle granular systems which are in contact with an external 
pressure bath. Our main concern was to discuss how constraint conditions can 
be applied to determine the interaction between the particles in an Andersen- 
type of dynamics. We have presented the results of some numerical tests and 
discussed the effect of the main parameters on the efficiency of the simulations 
and on the physical properties of the final packings. 

We restricted our study to the simple case where we allow only spherical strain 
of the system in order to achieve the desired pressure. However, the method 
can be generalized to apply other type of constraints to the stress tensor and, 
consequently, to allow more general strain deformations where shape as well 
as size of the simulation cell can be varied [24,25]. 
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Fig. 1. Two rigid particles before a possible contact. 




Fig. 2. The flowchart of the force calculation of a single contact in contact dynamics. 
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Fig. 3. Typical time evolution of the dilation rate A (left) the inner pressure P\ n 
(middle) and the system size L (right) during the compression of a 2D polydisperse 
sample. The data shown here were recorded in the reference system specified in the 
text. 
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Fig. 4. CPU time versus the simulation parameters. 
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Fig. 5. Mean overlap in terms of the simulation parameters. 
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Fig. 6. CPU time in terms of the mean overlap. The different curves are obtained by 
the variation of different parameters according to Figs. 4 and 5. These parameters 
are the external pressure P e xt ; the number of iterations per time step Nj, the length 
of the time step At (left) and the inertia of the system M\ (right). The open circle 
on the right indicates the most efficient simulation we could achieve by controlling 
the inertia M\ . 
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Fig. 7. Volume fraction versus the simulation parameters. 
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Fig. 8. (left) The angular distribution of the contacts. The number of contacts is 
plotted in terms of the direction of their normal vectors, (right) Spatial distribution 
of the contacts. The system contains N = 1000 disks in both figures. 
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